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Abstract. We present here a new solution for the astronomical computation of the orbital motion of the Earth 
spanning from to —250 Myr. The main improvement with respect to the previous numerical solution La2004 

(Laskar ot al. 2004) is an improved adjustment of the parameters and initial conditions through a fit over 1 Myr 
to a special version of the high accurate numerical ephemeris INPOP08 (Fienga et al. 2009). The precession 
equations have also been entirely revised and are no longer averaged over the orbital motion of the Earth and 
Moon. This new orbital solution is now valid over more than 50 Myr in the past or in the future with proper phases 
of the eccentricity variations. Due to chaotic behavior, the precision of the solution decreases rapidly beyond this 
time span, and we discuss the behavior of various solutions beyond 50 Myr. For paleoclimate calibrations, we 
provide several different solutions that are all compatible with the most precise planetary ephemeris. We have 
thus reached the time where geological data are now required to discriminate among planetary orbital solutions 
beyond 50 Myr. 



1. Introduction 

Due to gravitational planetary perturbations, the elliptical 
elements of the orbit of the Earth are slowly changing in 
time, as is the orientation of the planet's spin axis. As 
described by (Milankovitch 1941) these changes induce 
variations of the insolation received on the Earth's sur- 
face that are at the origin of large climatic changes. Since 
the work of (Hays et al. 1976), that established a corre- 
lation between astronomical forcing and the d^^O records 
over the past 500 kyr, there has been a increasing need 
for precise long term ephemeris for the Earth orbital and 
rotational evolution (see Laskar et al. (2004) for a more 
detailed historical account). 

For paleoclimate studies, the most widely used orbital 
solutions are nowadays either the averaged solution of 
(Laskar 1988; Laskar et al. 1993b) or the most recent nu- 
merical solution of Laskar et al. (2004). 

The first long term direct numerical integration (with- 
out averaging) of a realistic model of the Solar system, 
together with the precession and obliquity equations, was 
made by (Quinn et al. 1991) over 3 Myr. Over its range, 
this solution presented small differences with the secular 
solution of (Laskar 1988, 1990), (see Laskar et al. (1992)). 
The orbital motion of the full Solar system has then been 
computed over 100 Myr by (Sussman & Wisdom 1992), us- 
ing a symplectic integrator with mixed variables (Wisdom 
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& Holman 1991), confirming the chaotic behavior found by 
Laskar (1989, 1990). Following the improvement of com- 
puter technology, long term integrations of realistic mod- 
els of the Solar system have improved (Varadi et al. 2003; 
Laskar et al. 2004), but the main limitation remains the 
exponential divergence of nearby orbits resulting from the 
chaotic motion of the Solar system (Laskar 1989, 1990, 
1999). Although it is now possible to integrate the mo- 
tion of the Solar system over time periods of more than 5 
Gyr, comparable to its age or expected life time (Laskar 
& Gastineau 2009), it is clear that the chaotic behavior 
of the solution will still limit its validity to a few tens of 
Myr. 

The present paper is a contimiation of the work that 
has been conducted for decades in our group in order to 
obtain the most precise solution for the past evolution 
of the orbit and rotational state of the Earth, aimed to 
paleoclimate studies. 

The numerical integrator is the same symplectic inte- 
grator from (Laskar & Robutcl 2001) as the one used in the 
La2004 solution (Laskar et al. 2004), but it was entirely 
rewritten in C in order to access to extended precision on 
the Intel architecture. On the other hand, the tidal model 
has been largely modified, and is now close to the one used 
in the JPL planetary ephemeris DE405 (Standish 1998b) 
or in our new planetary ephemeris INPOP (Fienga et al. 
2008, 2009). The precession equations for the evolution of 
the spin axis of the Earth are also new (Boue & Laskar 
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2006). There are no longer averaged over the orbital mo- 
tion of the planets, and allow a precise computation of the 
evolution of the Earth spin axis that can be compared to 
the most precise model adopted by the lAU (Soffel et al. 

2003) (see Fienga et al. 2008). 

In previous long term solutions Laskar (1988); Quinn 
et al. (1991); Laskar et al. (1993a); Varadi et al. (2003); 
Laskar et al. (2004), the initial conditions of the solutions 
were obtained either directly from a high precision planet 
ephemeris, or by a fit over its full time span (as in La2004) 
that was still limited to a few thousands of years. It was 
thus difficult to monitor the real uncertainty of the used 
ephemeris. 

In the present work, we have profoundly removed this 
limitation. Indeed, in the past few years, we have entirely 
build a new high precision planetary ephemeris INPOP 
that has been fitted to all available planetary and Lunar 
observations (Fienga et al. 2008, 2009). This ephemeris 
that has been released already in two versions (INPOP06 
and INPOP08) is thus equivalent to the JPL ephemerides 
DE that were used for determining the initial conditions 
of the previous long term solutions. In addition, we have 
removed in INPOP all time limitations and carefully de- 
signed the numerical integrator. We could thus extend the 
integration of INPOP06 and INPO08 over 1 Myr with the 
full ephemeris model. The initial conditions of the present 
long term ephemeris could then be fitted over an extended 
interval of several hundreds of thousands of years before 
being extended to 250 Myr. By doing so, we were able 
to take into account the full precision of the ephemeris, 
and it appears now that the limitations is no longer in the 
model but in the planetary observations themselves. 

The first sections of this paper (sec 2 to sec 6 ) de- 
scribe successively the La2010 numerical model, its link 
with the INPOP ephemeris, the various La2010 solutions 
and their comparison with the high precision INPOP 
ephemeris and the previous La2004 solution. The follow- 
ing sections (7 to 9) are focusscd on the long term cy- 
cles that are present in the eccentricity solution and their 
stability. This topic is of essential importance for the at- 
tempt to establish an astronomically calibrated geological 
timescale (see Piilikc & Hilgen 2008) . Indeed, the La2004 
solution (Laskar et al. 2004) has been successfully used 
for the astronomical calibration of the Neogene period 
(w 23 Myr) (Lourens et al. 2004) that is included in the 
most recent standard timescale GTS2004 (Gradstein et al. 

2004) . At present, there is a continuous effort to improve 
this timescale and to extend the astronomical calibra- 
tion to the full Cenozoic period (« 65 Myr) through the 
Earthtime and Earthtime-eu projects (http://www. earth- 
time. org, http://earthtime-eu.eu). In order to do so, the 
length of validity of the orbital solution has to be extended 
by more than 20 Myr. Because of the chaotic behavior of 
the solution (Laskar 1989), this corresponds to improve 
the precision of the model and parameters by two orders 
of magnitude, and the present work is an attempt in this 
direction. 



On the other hand, beyond the horizon of predictibility 
of the orbital solution, it is tempting to use the recorded 
geological information to provide constraints on the or- 
bital motion of the Solar System (Lourens et al. 2001; 
Palike et al. 2004). Due to the lack of precision of the ge- 
ological data in remote periods of time, this can only be 
done through macroscopic aspects of the orbital solution. 
The analysis of the secular resonance 54 — 53 — 2(s4 — S3) 
(Sec. 8) is devoted to this problem. In particular, it is 
shown how the analysis of the modulation of the amplitude 
of the 405 kyr eccentricity term can discriminate among 
various orbital solutions and thus provide feedback from 
geological data to astronomical models. 

2. Numerical model 

The orbital solutions La90-93 (Laskar 1990; Laskar et al. 
1993a) were obtained by a numerical integration of the 
averaged equations of the Solar system, including the main 
general relativity and Lunar perturbations. As computer 
technology now allows to integrate directly precise models 
of the evolution of the Solar System over several hundreds 
of Myr, we have decided, since La2004 (Laskar et al. 2004) 
to use direct integrations, without averaging. 

The dynamical model and numerical integrators are 
very close to the ones of La2004. We will thus refer to 
(Laskar et al. 2004) for a detailed description of these 
models, and will only report here the elements that are 
different in the present model and integration. 

2.1. Dynamical model 

The orbital model comprises all 8 planets of the Solar 
system, and Pluto. The post-Newtonian general relativity 
corrections of order 1/c^ due to the Sun are included fol- 
lowing Saha & Tremainc (1994). The Moon is treated as 
a separate object. In order to obtain a realistic evolution 
of the Earth-Moon system, we also take into account the 
most important coefficient in the gravitational potential 
of the Earth and of the Moon, and the tidal dissipation in 
the Earth-Moon System (Laskar et al. 2004). Contrary to 
La2004, the precession and obliquity are now integrated 
without averaging over the orbital periods following (Boue 
& Laskar 2006). In the final runs, we have also added the 
contribution of the main 5 minor planets and some small 
correction to precession motions are made to take into 
account the remaining asteroids or other unmodelized pa- 
rameters. 

2.2. Numerical integrator 

As in La2004, the numerical integration was performed 
with the symplectic integrator scheme SABAC4 of 
(Laskar & Robutel 2001), with a correction step for the 
integration of the Moon. This integrator is particularly 
adapted to perturbed systems where the Hamiltonian gov- 
erning the equations of motion can be written as the sum 
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of an integrable part (the Keplerian equations of the plan- 
ets orbiting the Sun and of the Moon around the Earth), 
and a small perturbing potential representing the interac- 
tions among the planets. 

The step size used in the integration is r = 5 x 
10-3yr = 1.82625 days, while for La2010a, t = IQ-^yr = 
0.36525 days. The initial conditions of the integration were 
least square adjusted to a special version of INPOP that 
has been extended in time over 1 Myr. Depending on the 
solution, this fit was performed over 1 Myr or 580 kyr (see 
Tab. 2). 

In La2004, the integration was made in double pre- 
cision, with machine Epsilon £m ~ 2.22 x 10~^^. Here, 
we have integrated the solutions in extended precision on 
Xcon Intel processors, which allows arithmetics in 80 bits 
instead of 64 bits in double precision. The machine Epsilon 
becomes then e'^ « 1.1 x 10~^^. 

The integration time for our complete model, includ- 
ing 5 asteroids and the Moon as a separate object with 
T = 5 X 10~^yr is about one day per 3 Myr in extended 
precision, and one day per 6 Myr in double precision on 
a Intel Xeon E5462 2.8 Ghz workstation. When the step 
size is decreased to r = 10~^yr, for the nominal solution 
La2010a, the requested time is 5 days for 3 Myr, and more 
than one year for the whole integration. 

2.3. Numerical error 

As in La2004, the numerical error is estimated by com- 
paring two integrations with the same model and slightly 
different step size. For the nominal solution, we use 
T = 5 X lO""^ yr, and for the alternate solution r* = 
4.8828125 x 10~^ years. This special value is chosen in 
order that our output time span h = 1000 years cor- 
responds to an integer number (204800) of steps, in or- 
der to avoid any interpolation problems in the check of 
the numerical accuracy. With t = 10~^yr, we have then 
T* = 0.9765625 x lO-^yr (Tab. 2). 

2.3.1. Rotational evolution 

Contrarily to La2004, the precession equations are no 
longer averaged over the orbital motion of the planets 
or the Moon, but are treated in a vectorial manner, fol- 
lowing (Boue & Laskar 2006). Following (Darwin 1880; 
Mignard 1979), we assume that the torque resulting from 
tidal friction is proportional to the time lag At needed for 
the deformation to reach the equilibrium. This time lag is 
supposed to be constant, and the angle between the direc- 
tion of the tide-raising body and the direction of the high 
tide (which is carried out of the former by the rotation of 
the Earth) is proportional to the speed of rotation. Such 
a model is called "viscous", and corresponds to the case 
for which 1/Q is proportional to the tidal frequency. 

Various additional small dissipative effects as core- 
mantle friction (Poincare 1910; Rochester 1976; Lumb & 
Aldridge 1991; Correia et al. 2003), atmospheric tides 



(Chapman & Lindzen 1970; Volland 1978; Correia & 
Laskar 2003), mantle convection (Forte & Mitrovica 1997), 
climate friction (Rubincam 1990, 1995; Bills 1994; Ito 
et al. 1995; Levrard & Laskar 2003), have been discussed 
in La2004, but their effects are considered to be too small 
and too uncertain to be added in the model, as it was the 
case for La2004. 

3. The numerical ephemeris INPOP 

The initial conditions of La2004 were obtained by adjust- 
ment to the JPL numerical ephemeris DE406 (Standish 
1998a) over the full range of DE406, that is from -5000 
yr to -1-1000 yr from the present date. DE406 is itself ad- 
justed to planetary observations. 

With this procedure, we are limited by the range of the 
available ephemeris, and in general, the latest ephemeris 
is not always computed over a long time interval. For 
example, the most recent ephemeris from JPL, DE421 
(Folkner 2008), has only been provided over the time in- 
terval [1900, 2050] yr. Moreover, it is difficult to estimate 
the true uncertainty of the provided ephemeris. Most of- 
ten, this uncertainty is only revealed with the publication 
of the new ephemeris that can be the compared with the 
previous one. In order to overcome these limitations, we 
have undertaken in our group the construction of a full size 
planetary and lunar ephemeris. After five years of work, 
the solutions are now mature and two successive versions 
have already been published : INPOP06 (Fienga et al. 
2008) and INPOP08 (Fienga et al. 2009). The detailed in- 
formation about the dynamical models and fit to available 
observations can be found in the related publications. 

We have removed in the construction of these 
ephemerides all elements that would limit the length of 
validity of the solutions. In particular, we have not used 
some precession formulas for the evolution of the spin axis 
of the Earth. Instead, we have integrated together with 
the full ephemeris, a precession model for the Earth that 
is obtained after averaging over the rotation period of the 
Earth, but not over the orbital period of the Earth or of 
the Moon (Fienga et al. 2008). 

The full ephemeris could then be prolongated over 1 
Myr using extended precision 80 bits arithmetics with the 
Adams integrator of INPOP. This integration took about 
4 month of CPU on an Itanium 9040 1.6 Ghz workstation. 
This process was first made for INPOP06, and then for 
INPOP08, when the final version of this latest ephemeris 
(Fienga et al. 2009) was finally made available. These 
highly accurate ephemerides are then used for the cali- 
bration and evaluation of the long models La2010. 

We refer to (Fienga et al. 2008, 2009) for a precise 
description of INPOP06 and INPOP08. With respect to 
INPOP06, INPOP08 benefitted from several additional 
sets of observations. The Mars Express and Venus Express 
ranging data provided very precise measures of Earth- 
Mars and Earth- Venus distances with a precision of a few 
meters (Fienga et al. 2009). For Mars, this was a contin- 
uation of a long sequence of very precise measures that 
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Table 1. Maximum diflcrcnce between INPOP06 and 
INPOP08 over 1 Myr (top), and between INPOP08 (computed 
in extended precision) and INPOPOSd, a version of INPOP08 
computed in double precision, over 1 Myr (bottom). 





INPOP06-INPOP08 






a X 10" 


e X 10" 


i X 10" 


Mercury 




0.027 


25.544 


12.598 


Venus 




0.543 


4.746 


6.299 


EMB 




1.067 


4.709 


4.029 


Mars 




3.852 


7.134 


2.570 


Jupiter 




56.126 


20.542 


0.577 


Saturn 




585.092 


76.138 


1.315 


Uranus 




885.497 


92.313 


1.481 


Neptune 




3449.727 


104.593 


0.882 


Pluto 




19900.821 


297.483 


46.579 


Moon 




57.345 


42006.038 


1463.465 




INPOP08-INPOP08d 






a X 10" 


e X 10" 


i X 10" 


Mercury 




0.012 


0.077 


0.011 


Venus 




0.640 


1.706 


0.138 


EMB 




0.889 


1.768 


0.125 


Mars 




6.899 


6.081 


0.394 


Jupiter 




1.004 


0.175 


0.009 


Saturn 




2.542 


0.243 


0.007 


Uranus 




6.834 


0.333 


0.006 


Neptune 




13.479 


0.415 


0.008 


Pluto 




25.332 


0.455 


0.061 


Moon 




29.594 


12400.799 


527.701 



Note. EMB is the Earth-Moon barycenter. The semi-major 
axis a is in AU and the inclination with respect to the invariable 
plane i is in radians. 



had been acquired with the Martian spacecrafts since the 
first Viking landers on Mars, but for Venus, the new rang- 
ing data processed by ESOC were the first highly accu- 
rate estimates of the Earth- Venus distance which uncer- 
tainty was thus reduced from a few hundred meters to 
a few meters (Fienga et al. 2009). Another improvement 
in INPOP08 consists in the use of some Cassini normal 
points (Folkner 2008) that also help to constrain the po- 
sition of Saturn. In addition, in INPOP08, the Lunar or- 
bit was fitted to Lunar laser ranging data in a consistent 
way, while in INPOP06, the fit of the Lunar ephemeris 
was only made with respect to Lunar distances given by 
DE405 (Standish 1998a). 

It is always difficult to estimate the true uncertainty 
of an ephemeris. In (Fienga et al. 2009), this estimate is 
obtained by comparison with INPOP06 and DE421 over 
10 and 100 years. The differences between INPOP08 and 
DE421 are in general smaller, but comparable with the dif- 
ferences INPOP08-INPOP06, with the notable exception 
of the positions of Saturn, where the differences INPOP08- 
DE421 are one order of magnitude smaller than the dif- 
ferences INPOP08-INPOP06. This is certainly the conse- 
quence of the use in both DE421 and INPOP08 of the 
new Cassini data that constrain very much the position of 
Saturn. 



After 100 yr, the differences INPOP08-DE421 in 
barycentric positions range from a few kilometers for the 
inner planets to a few thousands km for the outer plan- 
ets, and only 40 km for Saturn (Fienga et al. 2009, Tab. 
6). This is several order of magnitude more than the error 
coming from the numerical integration (Fienga et al. 2008, 
Tab. 1) that reach only a few micrometers after the same 
range, and less than a few meters after 10000 yr. It can 
thus be assumed that after one million year, the numeri- 
cal error in the integration of INPOP will still be smaller 
than the propagation of the uncertainty of the model and 
parameters, obtained by the fit to planetary positions. 

We have thus prolongated the two INPOP ephemeris 
(INPOP06 and INPOP08) over 1 Myr in order to use these 
solutions as a starting point for the long term ephemeris. 
The accuracy of these solutions after 1 Myr is then eval- 
uated by the comparison of INPOP06 to INPOP08, with 
the assumption that the real uncertainty of INPOP08 will 
be smaller than the difference INPOP08-INPOP06 (Table 
1). In Table 1, we also provide the differences of two in- 
tegrations of INPOP08 made in doTible prc^cision and in 
extended precision, in order to evaluate the numerical pre- 
cision of the integration. Prom the comparisons made in 
(Fienga et al. 2008), we can assume that the error on the 
integration of INPOP08 made in extended precision is in 
fact several orders of magnitude smaller than the differ- 
ences reported in this table. 

4. Successive versions of La2010 

The process leading to a long term solution is long, as we 
had to wait first for the INPOP solution to be ready over 1 
Myr, and then only we could make the fit of the long term 
model. After that, the integration of the long term model 
over 250 Myr in extended precision still required about 
3 months of CPU time, and more than one year when 
the step size is reduced to 10^'^ yr. This is why we have 
performed several versions of these long term ephemeris, 
that could be used for comparisons, and also to study the 
stability of the solution with respect to improvement of 
the INPOP ephemeris. As the INPOP08 ephemeris was 
only finished very recently, some of the solutions fitted 
to INPOP08 have been fitted over only 580 kyr instead 
of 1 Myr for INPOP06, but this did not make a large 
difference, and the solutions are still at the end compared 
to INPOP08 over IMyr as the integration of INPOP08 
has now reached IMyr. The various models that have been 
selected are summarized in Table 2. 

The solution La2010d has been fitted to INPOP06 over 
1 Myr, while the more recent solutions La2010a, La2010b, 
La2010c, and their associated solutions La2010a*, 
La2010b*, La2010c*, have been fitted over INPOP08, over 
580 kyr for the solutions with index a,b, and over 1 Myr for 
La2010c and La2010c*. All these models, except La2010c 
and La2010c* comprise the five major asteroids, Ceres, 
Vesta, Pallas, Iris and Bamberga. 

In all cases, the parameters are taken from the corre- 
sponding INPOP ephemeris, as well as the starting value 
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Table 2. Variants of the La2010 solutions. The nominal solution La2010a is obtained for a stepsize t = 1 x 10 in extended 
precision. The other solutions, with different setting have been computed to test the stability of the nominal solution. 



name 


files 


ephem 


fit 




prcc 




La2010a 


ast5AL08cxc 


INPOPOSa 


0.5S Myr 


1 X 10^^ 


Ext 


5ast,M 


La2010a* 


ast5AL08czc 


INPOPOSa 


0.58 Myr 


0.9765625 x 10^^ 


Ext 


5ast,M 


La2010b 


ast5AL08cx 


INPOPOSa 


0.58 Myr 


5 X IQ-^ 


Ext 


5ast,M 


La2010b* 


astSALOScz 


INPOPOSa 


0.58 Myr 


4.8828125 x 10"^ 


Ext 


5ast,M 


La2010c 


astOAL08cx2a 


INPOPOSa 


1 Myr 


5 X IQ-^ 


Ext 


Oast,M 


La2010c* 


astOAL08cz2a 


INPOPOSa 


1 Myr 


4.8828125 x IQ-^^ 


Ext 


Oast,M 


La2010d 


astSALix 


INPOP06 


1 Myr 


5 X 10"^ 


Ext 


5ast,M 



Note. The column "files" denotes the name of the computer files of the solution, "ephem" is the name of the reference INPOP 
ephemeris. "fit" is the time length of the ephemeris used for the fit. t is the step size of the numerical integration. In "prec", 
Ext indicates that the integration was performed in extended precision. In the last column, 5ast indicates that 5 asteroids have 
been integrated, and M stands for the Moon as a separate object. 




-100 -90 



-70 -60 -50 -40 -30 -20 -10 
time (Myr) 



Fig. 1. Estimate of the numerical precision of the solutions 
La2010a,b,c. The estimate is obtained by the difference in the 

eccentricity of the Earth obtained with the integration of two 
solutions La2010x and La2010x* for x = a,b,c (see Tab. 2). 



of the initial conditions. In order to take into account the 
differences of models, we then perform a fit of the semi 
major axis, and add a small precessing term that can be 
thought as representative of the average contribution of 
the minor planets that have not been taken into account 
in our simplified models. In these solutions, the Moon is 
integrated as a separate object, taking into account the 
tidal dissipation in the Earth Moon system. The step size 
is then 5 x 10~^ yr for La2010b,c,d and 10~^ yr for the 
nominal solution La2010a. 

In order to check the numerical accuracy of the solu- 
tion, we have also integrated these solutions with an al- 
ternate step size of r* = 4.8828125 x 10~^ yr for b and c, 



and T* = 0.9765625 x lO'^ yr forLa2010a*. These values 
are different, but close to the nominal step size. They are 
taken in such way that 1024 x r* = t, where t is the nom- 
inal step size of the corresponding solution. In Fig.l, the 
difference in the eccentricity of the Earth for two solutions 
La2010x and La2010x* arc plottcxl over tiuK^ for 100 Myr 
for X = a,b, c. Because of the exponential divergence of 
the solutions resulting from chaotic behavior, the differ- 
ence is nearly zero for a very long time;, of more than 50 
Myr, and then grows rapidly to maximal value, as the two 
solutions will become out of phase. 

This is an external way to evaluate the precision of 
the numerical integration, but it is in fact a pessimistic 
view. Indeed, we have fitted the solutions La2010a,b,c to 
INPOP08, but the initial conditions of La2010x* is the 
same as for La2010x. The difference of step size will then 
induce a difference of reference Hamiltonian in the sym- 
plcctic integration of the system, which should explain 
most of the difference that is observed here. This is why 
for a reduced step size, as in La2010a, the difference be- 
tween La2010a and La2010a* is smaller. 

Nevertheless, although pessimistic, this shows that 
the numerical error can be neglected over 55 Myr for 
La2010b,c and 60 Myr for La2010a. This is why we have 
selected La2010a as our nominal solution. 

5. Comparison with INPOP08 

The solution La2004 was fitted to DE406 over its full 
range, that is over the interval [—5000 : -1-1000] yr from 
now. In 2004, there was no possibility to compare it to an 
accurate ephemeris over a longer time. With the construc- 
tion of the new INPOP ephemerides, this becomes now 
possible, as we have extended INPOP06 and INPOP08 
over 1 Myr. As the set of observation used in INPOP08 
is significantly larger than the one of INPOP06, we will 
use INPOP08 as the reference ephemeris, representing the 
best knowledge of the orbital motion of the Solar System 
that we can achieve at present. 

We have thus compared La2004 to INPOP08 as well 
as INPOP06 and the new computed solutions of the Earth 
eccentricity for La2010a,b,c,d (Fig. 2). From Fig. 2 (top), 
it is clear that the new solution La2010a is a significant 
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-1.5e-06 I 1 1 1 1 1 

-1 -0.8 -0.6 -0.4 -0.2 

time (Myr) 



Fig. 2. Differences in the eccentricity of the Earth Moon 
barycenter over 1 Myr for various solutions as follows : 
al = La2004-INPOP08; a2 = La2010d-INPOP08; a3 =b3= 
La2010a - INPOP08; a4=INPOP06-INPOP08; bl = La2010b 
- INPOP08; b2 = La2010c -INPOP06. It should be noted that 
the vertical scale is enlarged 10 times in the bottom plot. All so- 
lutions are compared to INPOP08. On the top figure, La2010d 
(a2) and INPOP06 (a4) are almost superposed. This is because 
La2010d was adjusted over INPOP06. 



improvement with respect to La2004. Indeed, the differ- 
ence La2010a-INPOP08 (Fig. 2 (a3) and (b3)) is nearly 15 
times smaller than the difference La2004-INPOP08 (Fig. 
2 (al)). 

In Fig. 2, we can see that La2010d-INPOP08 (a2) is 
almost superposed with INPOP06-INPOP08 (a4). This is 
because La2010d has been adjusted to INPOP06. It is 
also clear from this plot that the differences between a 
long time solution and its reference ephemeris are now 
much smaller than those of two consecutive versions of 
the high resolution planetary ephemeris (as INPOP06 and 
INPOP08). 

The main result of these comparisons, are also dis- 
played in Tables 3 and 4, for the various solutions that we 
have selected. The differences between the long term "sim- 
plified" model that we use here and the most precise plane- 
tary ephemerides are now much smaller than the difference 





ast5AL08cx -INPOP08 1 Ma 




a X 10" 


e X 10" 


i X 10" 


Mercury 


0.006 


0.583 


0.498 


Venus 


0.341 


1.027 


0.351 


EMB 


0.453 


1.182 


0.343 


Mars 


1.774 


1.608 


0.351 


Jupiter 


0.551 


0.147 


0.049 


Saturn 


1.758 


0.359 


0.037 


Uranus 


5.129 


0.478 


0.029 


Neptune 


8.734 


0.306 


0.027 


Pluto 


17.901 


0.325 


0.049 


Moon 


8.916 


22005.375 


12538.037 




ast5SL08ax-INPOP08 1 Ma 




a X 10'' 


e X 10" 


i X 10" 


Mercury 


0.002 


3.962 


2.859 


Venus 


0.155 


1.975 


0.978 


EMB 


0.224 


1.540 


0.968 


Mars 


1.744 


2.823 


0.954 


Jupiter 


0.278 


0.113 


0.033 


Saturn 


0.777 


0.273 


0.023 


Uranus 


2.256 


0.352 


0.020 


Neptune 


4.052 


0.147 


0.009 


Pluto 


8.805 


0.165 


0.019 



Table 3. Maximum difference between the La2010 solu- 
tions and INPOP08 over 1 Myr. Top: Maximum difference 
ast5AL08cx -INPOP08 over 1 Myr. Bottom : Maximum dif- 
ference ast5SL08ax -INPOP08 over 1 Myr. In ast5SL08ax, the 
Moon contribution is averaged. All solutions are in extended 
precision. EMB is the Earth-Moon barycenter. The semi-major 
axis a is in AU and the inclination with respect to the invari- 
able plane i is in radians. 



between two consecutive planetary ephemeris (INPOP06- 
INPOP08) that can be considered as representative of the 
true uncertainty of the ephemeris. The main limitation of 
the precision of the long term planetary solution then re- 
sides in the precision of the planetary ephemeris, that is 
in the planetary observations. 

6. Comparison with La2004 

After comparing the solutions over 1 Myr with the most 
precise planetary ephemeris, we will now compare the vari- 
ous solutions La2010a,b,c,d to the former solution La2004 
over the whole expected range of validity of these solu- 
tions, that is over a few tens of million of years. 

In (Fig. 3), we have represented the variation of the ec- 
centricity of the EMB from -30 Myr to -50 Myr for the 
La2004 solution and the four La2010a,b,c,d new solutions. 
The interval [—30 : 0] Myr is not represented as all five so- 
lutions practically coincide over this time interval. Indeed, 
some discrepancies between La2004 and the new La2010 
solutions only appear on the time interval [—40 : —30], 
although most of the time the solutions are still very sim- 
ilar, and the small differences that can be seen on Fig. 3 
will most probably not lead to any significant change in 
the paleoclimate records. We can even consider that the 
solution La20G4 is still in good agreement with the new 
solutions until —45 Myr. This is in good agreement with 
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Fig. 3. Eccentricity of the Earth for the solutions La2004 (L04), La2010a (LlOa), La2010b (LlOb), La2010c (LlOc), and La2010d 
(LlOd). On the interval [—40 : —30] Myr, the four solutions are virtually identical, but before —45 Myr, La2004 begins to depart 
significantly from the La2010 solutions. We have then plotted again the time interval [—50 ; —40] Myr (bottom), removing 
La2004 for a better comparison of the La2010 solutions. Over this interval, these three solutions are very similar. 



the expected precision that was forecasted in (Laskar et al. 
2004). 

Beyond —45 Myr, noticeable differences become to ap- 
pear, and the solution La2004 becomes significantly dif- 
ferent than the La2010 solutions. We have thus made an 
additional plot of the [—50 : —40] Myr interval on Fig. 3, 
with only the solutions La2010a,b,c,d. These latest solu- 
tions well agree on this time interval, despite the variations 
of models or initial conditions among these new solutions. 



We can thus consider that the present new solutions are 
at least valid over 50 Myr. 

Beyond —50 Myr, the situation is more confused as the 
solutions La2010a,b,c,d present significant variations (Fig. 
4). Nevertheless, from Fig. 4, it can be seen that if moder- 
ate precision is only required, the solution could eventually 
used up to —60 Myr. In particular, the solutions are still 
well in phase around —55 Myr. This date is of particular 
interest, as it corresponds to specific climatic events that 
have been well documented in the paleoclimate records : 
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-63 -62 
time (Myr) 



Fig. 4. Eccentricity of the Earth for the solutions La2010a (LlOa), La2010b (LlOb), La2010c (LlOc), and La2010d (LlOd) from 
—50 to —65 Myr. Although the various solution begin to diverge beyond 53 Ma, it is remarkable that the minima of eccentricity 
at 51.75 Ma, 56.25 Ma and 58.25 Ma (in grey) correspond in all various solutions. 



the Paleocene-Eocene Thermal Maximum (PETM) that 
is dated aromid 55.53-56.33 Ma^ (Westcrhold ct al. 2007, 
2008), and the Elmo Thermal Maximum (ETM) dated at 
53.7-54.5 Ma (Lourens et al. 2005; Westerhold et al. 2007, 
2008). 

On the other hand, from Fig. 4, it is quite clear that 
the solutions La2010a,b,c,d cannot be used beyond 60 



^ Myr design the duration of 1 million of years, while Ma 
stands for mega-annum, and represent date in the past from 
the present. 



Ma, as the solutions have large differences in the inter- 
val [—65 : —60] Myr. It should thus be stressed that in 
this time interval, the direct use of the eccentricity solu- 
tion of the Earth for geological calibration should be used 
with utmost care. 



Practically, if one wishes to use these solutions be- 
yond 50 Myr, for a geological calibration for example, one 
should use only the features of the solutions that remain 
the same in the four La2010 solutions. This will assess the 
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astSALix -INPOP06 1 Ma 




a X 10" 


e X 10" 


i X 10" 


Mercury 


0.005 


0.473 


0.342 


Venus 


0.034 


0.503 


0.209 


EMB 


0.049 


0.505 


0.234 


Mars 


0.721 


0.873 


0.397 


Jupiter 


0.121 


0.069 


0.046 


Saturn 


0.980 


0.157 


0.032 


Uranus 


1.553 


0.139 


0.025 


Neptune 


2.263 


0.068 


0.020 


Pluto 


7.068 


0.104 


0.022 


Moon 


21.439 


26686.705 


12769.784 




ast5ALh-INPOP08 1 Ma 




a X 10" 


e X 10" 


i X 10" 


Mercury 


0.007 


0.503 


0.357 


Venus 


0.290 


1.156 


0.225 


EMB 


0.379 


0.986 


0.257 


Mars 


2.757 


2.484 


0.469 


Jupiter 


0.563 


0.139 


0.047 


Saturn 


2.406 


0.272 


0.034 


Uranus 


5.415 


0.268 


0.026 


Neptune 


9.845 


0.320 


0.021 


Pluto 


19.768 


0.361 


0.052 


Moon 


21.093 


26481.398 


12727.162 



Table 4. Maximum difference between the La2010 solutions 
and INPOP06 over 1 Myr. Top: Maximum difference ast5ALix 
-INPOP06 (Both solutions are in extended precision) over 1 
Myr. Bottom: Maximum difference ast5ALh -INPOP06 over 1 
Myr. The solution ast5ALh has been computed in double pre- 
cision, while INPOP06 is in extended precision. Both solutions 
have been fitted to INPOP06 over 1 Myr. EMB is the Earth- 
Moon barycenter. The semi-major axis a is in AU and the 
inclination with respect to the invaraible plane i is in radians. 



Table 5. Frequency decomposition of z = eexpiru for the 
Earth on the time interval [—15, 4-5] Myr (Laskar et al, 2004). 



n 




Mfc C'/yr) 


bk 


(fik (degree) 


1 


95 


4.257564 


0.018986 


30.739 


2 


92 


7.456665 


0.016354 


-157.801 


3 


9i 


17.910194 


0.013055 


140.577 


4 


93 


17.366595 


0.008849 


-55.885 


5 


51 


5.579378 


0.004248 


77.107 



stability of such a calibration with respect to the uncer- 
tainty of the La2010 solutions. 

7. Long term cycles 

The complete eccentricity solution of the Earth allow a 
direct adjustment of paleoclimate data to the oscillations 
of about 95 kyr and 124 kyr of the eccentricity (see Laskar 
et al. 2004), but for ancient records, this signal may not 
be clearly visible in the sediments. On the other way, 
the 405 kyr oscillation with argument 52 ~ 95, where 
(Table 6) are the secular frequencies of the Solar System 
(see Laskar et al. 2004) is very often present in the sedi- 
mentary records. This term is the largest term in a quasi 
periodic approximation of the eccentricity of the Earth 



solution of the Earth 9 




1 2 3 4 5 

freq "/yr 



Fig. 5. In (a) is the spectrum of the eccentricity over 65Myr 
from the La2010a solution, limited to the interval [0,5" /yr]. 
In (b), the same spectrum is plotted for a solution build with 
only the 5 main terms of 2:3 (Table 5). The main peaks are 
then easily identified and thus also the main peaks of the full 
eccentricity (a). 

Table 6. Main secular frequencies gi and Si of La2004 and 
La2010a determined over 20 Ma for the four inner planets, and 
over 50 Ma for the 5 outer planets (in arcsec yr^^). Aioo are 
the observed variations of the frequencies over 100 Myr (Laskar 
et al. 2004). In the last column, the period of the secular term 
are given. 





La2004 


La2010a 


Aioo 


period (yr) 


91 


5.59 


5.59 


0.13 


231 843 


92 


7.452 


7.453 


0.019 


173 913 


9 s 


17.368 


17.368 


0.20 


74 620 


9i 


17.916 


17.916 


0.20 


72 338 


95 


4.257452 


4.257482 


0.000030 


304 407 


56 


28.2450 


28.2449 


0.0010 


45 884 


97 


3.087951 


3.087946 


0.000034 


419 696 


58 


0.673021 


0.673019 


0.000015 


1925 646 


59 


-0.34994 


-0.35007 


0.00063 


3 703 492 


Sl 


-5.59 


-5.61 


0.15 


231 843 


S2 


-7.05 


-7.06 


0.19 


183 830 


S3 


-18.850 


-18.848 


0.066 


68 753 


S4 


-17.755 


-17.751 


0.064 


72 994 


S5 
Sd 


-26.347855 


-26.347841 


0.000076 


49188 


S7 


-2.9925259 


-2.9925258 


0.000025 


433 079 


S8 


-0.691736 


-0.691740 


0.000010 


1 873 547 


S9 


-0.34998 


-0.35000 


0.00051 


3 703 069 



(see Laskar et al. 2004, Table 6), and is less influenced by 
the chaotic diffusion present in the Solar system than the 
shorter period terms around 100 kyr (Laskar 1990; Laskar 
et al. 2004). 
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-26 -24 
time (Myr) 



Fig. 6. Filtered eccentricity of the Earth for the solutions La2004 (L04), La2010a (LlOa), La2010b (LlOb), La2010c (LlOc), and 
La2010d (LlOd) from +10 to -30 Myr. The solution is filtered in the interval [0 : 2.2"/yr], that is for periods in [589, +oo[ kyr. 



In recent works, the modulation of the 405 kyr compo- 
nent, due to the beat 53 — 174 of period ~ 2.4 Myr has also 
been identified in the sedimentary records, and is thought 
to be a key factor for the onset of special climate events 
(Lourens et al. 2005; Palike et al. 2006; van Dam et al. 



2006). More generally, there has been a large effort to 
search for long term cycles in the sedimentary records and 
to use them to hook the sedimentary data to the eccen- 
tricity computations (Olsen & Kent 1996; Lourens et al. 
2005; Westerhold et al. 2007, 2008; Jovane et al. 2010). 
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Fig. 7. Filtered eccentricity of the Earth for the solutions La2004 (L04), La2010a (LlOa), La2010b (LlOb), La2010c (LlOc), and 
La2010d (LlOd) from -30 to -70 Myr. The solution is filtered in the interval [0 : 2.2"/yr], that is for periods in [589, +oo[ kyr. 



In figure 5a is plotted the spectrum of the nominal so- 
lution La2010a, hmited to the range [0,5] "/yr (periods 
larger than 260 kyr). As there is a gap at about 2.2 "/yr 
in this spectrum, we have filtered the eccentricity data for 



all various solutions in the range [0,2.2] "/yr (Figures 6, 
7). 

Here again, it is clear that all solutions La2004, and 
La2010a,b,c,d are practically identical over [—30 : 0] Myr, 
and very similar up to - 45 Myr, where La2004 starts to 
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differ notably from the new solutions La2010, which still 
behave in a similar manner up to about -50 Myr where 
the situation becomes more confused. 



7.1. The modulation of the §2 — 405 kyr cycle 

In order to examine more closely the long term cycles in 
the Earth eccentricity, we have identified the origin of the 
main spectral terms in the eccentricity spectrum of fig- 
ure 5a. In order to do so, a synthetic eccentricity curve is 
build along the same time range using only the five terms 
of eexp(in7), as provided by the frequency decomposition 
of the solution La2004 over the time interval [—15,-1-5] 
Myr, taken from (Laskar et al. 2004). The plot of the spec- 
trum of the eccentricity function that is obtained with this 
purely quasiperiodic signal with frequencies limited to the 
linear terms 51,52,53,34,55 (Table 5) is plotted in Figure 
5b. 

As this synthetic model is quasiperiodic with only five 
main frequencies, the identification of the main spectral 
terms of the eccentricity are then obtained unambiguously 
with a spectral analysis over 65 Myr and are detailed in 
figure 5b. These terms are easily related to the correspond- 
ing peaks of the full eccentricity spectrum of Figure 5a. 

This exercise, that is in some sense complementary 
from a full quasiperiodic decomposition of the eccentric- 
ity as in (Laskar et al. 2004, Table 6) allows to better 
understand the behavior and origin of the main long term 
cycles observed by the practitioners when comparing to 
geological data (Olsen & Kent 1996; Lourens et al. 2005; 
Westerhold et al. 2007, 2008; Jovane et al. 2010; Hilgen 
et al. 2010). 

The leading periodic term is the well known 405 kyr 
term 52 — 55, but this term is surrounded by two terms 
(52 -55) -(54 -53) and (52 -55) + (54 -53) that wiU induce 
with 52 — 55 a modulation of the 405 kyr eccentricity term 
with a frequency of 54 — 53, corresponding to a 2.4 Myr 
period. An obvious consequence is that when analyzing 
geological data to search for the 405 kyr term, one needs to 
use a spectral window that includes these two side terms, 
that is a window similar to the [2.2, 4.3] " /yi window used 
in figure 8. Additionally, the two terms 51 — 55 of period 
« 1 Myr, and 52—51 of period « 688 kyr are also of strong 
amplitude in the eccentricity spectrum. 



7.2. The 54 - gs 2.4 Myr cycle 

Moreover, 54 — 53 appears also directly as a main periodic 
term of the eccentricity (Fig. 5a, b). The same 2.4 Myr cycle 
can thus be directly retrieved from the eccentricity curve. 
Indeed, in figure 8, we have plotted both the filtered ec- 
centricity in the interval [2.2,4.3] "/yr (e°) (in red) and 
as well the filtered eccentricity with a [0,0.1]"/yr window 
(e*") (in blue). It can then be seen that the envelope (in 
green) of e" is almost identical to the opposite of e'' (Fig. 
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Fig. 8. Filtered eccentricity around the 405 kyr period for 
La2010a. In red is e", the filtered eccentricity in the band 
[2.2, 4.3]" /yi- ([301,589] kyr period), while e^ the filtered ec- 
centricity in the band [0,1.1] "/yr ( period > 1.18 Myr) is 
plotted in blue. The opposite (in pink) of the maximum en- 
veloppe of e" (in green) nearly coincide with e''. 

As a consequence, the two components e'^ and of 
the eccentricity need to be added in order to really eval- 
uate the component of the 2.4 Myr term (Fig. 9). In the 



resulting 



curve, the variation of the maxima are 



then attenuated while the minima variations are increased 
to about 0.02. The variations of the minima are in phase 
with the 54 — 53 term (plotted in blue in Fig. 9). 

The large size of these variations makes it then un- 
derstandable that a signature of these variations could be 
recorded in the sedimentary paleoclimate signal. Indeed, 
although the global mean annual insolation on earth varies 
as and thus is not much influenced by eccentricity vari- 
ations, this is not the case for seasonal variations. Indeed, 
if one considers a black body with uniform temperature at 
distance d of a star, using Stefan's law for the emission of 
a black body, one finds that its surface temperature T is 
proportional to d~^^^. In this case, the difference ST be- 
tween perihelion and aphelion temperature will be given 

by 



6T 
~T 



1 2ae 

2 a 



(1) 



A change of 0.02 in the eccentricity corresponds thus in 
this simplify model to a change of about 0.02 x 300 = 6K 
in the difference between perihelion and aphelion temper- 
atures. 

Because of the increasing importance of this 2.4 Myr 
component in some of the analysis of sedimentary records. 
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Fig. 9. e°' + e'' (see Fig. 8) for the La2010a solution is plotted in 
red. Its minimum are in phase with e**, the filtered eccentricity 
in the band [0, 1.1] "/yr (in blue). 



we have added here a detailed comparison of the filtered 
solution in the [0, l.l]"/yr interval in figure 10 for the time 
intervals [—55, —40] and [—65, —50] Myr time intervals. It 
should be noted that it is not necessary to compare the 
various orbital solutions in the [—40, 0] Myr time interval 
as they are practically identical in this range. 

As in the previous discussion, we can see that all curves 
are very similar until —45 Myr, while La2004 differs sig- 
nificantly beyond —45 Myr. This is why this solution is 
no longer plotted on the bottom plot of figure 10, which 
is displayed on the [—65, —50] time interval. This range is 
particularly critical, as it corresponds both to the location 
of the PETM (at about -55 Myr) and of the K/P bound- 
ary (at about —66 Myr) (Lourens et al. 2005; Westerhold 
et al. 2007, 2008). The various solutions begin to differ sig- 
nificantly beyond —53 Myr, but it can be remarked that 
the two maxima at about —57.3 Myr and —59.3 Myr agree 
for all four La2010 solutions, although they largely differ 
around —55 Myr. One could thus use these three peaks in 
order to attempt to fit a geological time scale beyond -50 
Myr, in the [-60, -50] Myr interval. 

In fact, in the La2010a solution, the 54 — (73 argument 
has a period of about 27r/2.664 w 2.36 Myr in the interval 
[—45, 0] Myr, but beyond —45 Myr, this period changes 
due to chaotic diffusion (Fig. 11). As this occurs at the 
border of the validity range of the solution, it is still diffi- 
cult to be sure of the real behavior of the 54 — argument 
beyond —45 Myr, and it will be necessary to confront these 
data to geological records to confirm the behavior of the 
solar system eccentricity solution. 



La2004 
La2010a 
La2010b 
La2010c 
La2010d 




-60 -58 -56 
time (Myr) 

Fig. 10. Filtered eccentricity in the band [0, 1.1] " /yr ( period 
> 1.18 Myr) for La2004 and La2010a,b,c,d. 



8. Resonant angles 

8.1. Secular resonances 

The previous discussion demonstrates the importance of 
the behavior of the 34 — 33 argument in the macroscopic 
aspect of the variations of the Earth's eccentricity, and 
thus its possible relation with the past climate on Earth. 

Indeed, using the secular equations, Laskar (1990, 
1992) demonstrated that the chaotic behavior of the Solar 
System arise from multiple secular resonances in the inner 
Solar System, and in particular, from the critical argument 
associated to 



= (s4 - S3) - 2(54 - gs) 



(2) 



where 33, 34 are related to the precession of the perihelion 
of the Earth and Mars, S3, S4 are related to the precession 
of the node of the same planets. This argument is presently 
in a librational state, but can evolve in a rotational state, 
and even move to libration in a new resonance, namely 



(s4 - S3) - (514 - 33) = 



(3) 



The argument 9 as well as the other important reso- 
nant argument {a — (51— 55) — (si— S2)) that was identified 
by Laskar (1990, 1992) as the origin of the chaotic behavior 
of the inner planets are plotted in figure 12 for all solu- 
tions La2004, La2010a,b,c,d. In all cases, transition from 
libration to circulation appear around —50 Myr, leaving 
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time (Myr) 

Fig. 11. Argument (in radians) related to g4 — gs ~ 2.664 T 
versus T (in Myr) for La2010a (top). Argument (in radians) 
related to S4 - S3 - 2 x 2.664 T versus T (in Myr) for La2010a 
(bottom). 



some uncertainty to the behavior of the solution beyond 
this date. 



8.2. Searching for some geological evidence of chaos 

The transition from libration to circulation of the resonant 
argument related to 9 — (54 — S3) — 2{g4 — g^) is directly 
linked with the chaotic diffusion of the orbital trajectories. 
Searching in the geological record for the evidence of such 
a transition would thus be an observational confirmation 
of the past evolution of the Solar System. 

As it appears in the previous sections, it becomes more 
and more difficult to obtain by numerical computations 
only the date of the first transition from libration to cir- 
culation for this resonant argument (Fig. 12). The di- 
rect observation of the individual arguments related to 
93-, gi^ S3, S4 is certainly out of reach. On the other hand, 
the argument 9 corresponds to a 2 : 1 resonance between 
the two secular terms 54 — and S4 — S3, both terms be- 
ing present in the sedimentary records. We have discussed 
about the importance of the 34 — 53 beat in the eccentric- 
ity solution. In a similar way, S4 — S3 appears as a beat of 
about 1.2 million of years in the solution of obliquity, as 
the result of the beat between the p+s^ and P + S3 compo- 
nents of the obliquity, where p is the precession frequency 
of the axis (see Laskar et al. 2004, Fig. 7). 

With the occurrence of these beats, the detection of 
the resonant state in the geological data becomes pos- 
sible. Indeed, the modulation of 1.2 Myr of the obliquity 
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Fig. 12. Resonant arguments (in radians versus time in Myr) 
cr = (gi - gs) - (si - S2) (top) and 6 = (s4. - S3) - 2(34 - 
33 ) (middle and bottom) for the different solutions La2004, 
La2010a,b,c,d. In the bottom plot, is diplayed at a greater 
scale. One can see that the La2010 solutions are very close up 
to 50 Myr. 



appears clearly in the spectral analysis of the paleoclimate 
record from Ocean Drilling Program Site 926, (Zachos 
et al. 2001). Moreover, using the ODP legs 154 and 199, 
Palike et al. (2004) could find some evidence that the crit- 
ical argument of 9 did not show a transition to circulation 
at 25 Myr, as in La93, but remained in libration over 30 
Myr, as in the La2004 solution. 

Searching for a transition of the (S4 — S3) — 2(34 — 53) 
resonance to the (S4— S3) — (54 — (73) resonance, as displayed 
in figure 12) is difficult, as it requires to obtain both a good 
signal in eccentricity (or precession) and in obliquity. It 
may be more direct to search only for a modulation of 
the 34 — 33 (or S4 — S3) period, as it appears in figure 11. 
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Fig. 13. Argument (in radians) related to g4 — gs ~ 2.664 T 
versus T (in Myr) for various orbital solutions (top). Argument 
(in radians) related to S4 — S3 — 2 x 2.664 T versus T (in Myr) 
for various orbital solutions (bottom). 



Indeed, the change in La2010a at —45 Myr of the §4 — 
slope, from a period of about 2.4 Myr to a period of about 
2 Myr also reflects the same chaotic transition. 

In fact, this change can be directly seen in the eccen- 
tricity record (Figs 9, 10), as it will induce a change of the 
time interval from two maxima in the filtered eccentric- 
ity (Figs 10) from a 2.4 Myr period to a 2 Myr period. If 
the 405 kyr g2 — 55 signal is well present in the geological 
data, this becomes then a macroscopic feature that can be 
detectable. Indeed, a local time scale can be established 
using the 405 kyr signal, and the modulation period of 
this signal should be the — g^ term. The transition is 
then obtained as a transition from 6 periods per beat to 5 
periods per beats. Such geological data could then be able 
to discriminate among the various La2010 solutions (Fig. 
13). It becomes therefore even more important to search 
for good sedimentary sections where the 405 kyr signal is 
well determined. 

9. Stability of the 52 — 95 405 l<yr cycle 

As it was stressed above, the 52 — 55 405 kyr argument 
has a particular importance in long time geological cali- 
bration, as it is present in many sedimentary records and 
its good stability (Laskar 1990) can allow to use it as a ref- 
erence time scale This argument is indeed visible in many 
sedimentary records of the Early Mesozoic (Olsen & Kent 
1999, and references therein). As in (Laskar et al. 2004), 
we have tested the stability of this argument over the full 




-150 -100 
time (Myr) 

Fig. 14. 405 kyr term in eccentricity. Maximum difference (in 
radians) of the argument 6g^-g^{t) — dg^-g^{Qi) of (72 — 55 in 
all solutions La2004, La2010a,b,c,d with respect to the linear 
approximation 6a{t) — 3.200" f where t is in yr. 

period of our integrations, that is over 250 Myr by com- 
parison of its evolution on all retained La2010 solutions 
and La2004 (14). The present values of g2 — g^ do not 
differ significantly from the value of La2004 (Laskar et al. 
2004). The frequency 52 — 55 is thus kept to its La2004 
value 



= 3.2007yr 
which corresponds exactly to a period 
P405 = 405000 yr. 



(4) 



(5) 



As seen in Fig 14, the maximum deviation obtained by 
comparing all solutions is about 2tt over 250 Myr, which 
correspond to a full cycle of 405 kyr after 250 Myr, as was 
given already in (Laskar et al. 2004). 

10. Discussion and future work 

The new orbital solutions of the Earth that are pre- 
sented here can be used for paleoclimate computations 
over 50 Myr. Beyond that time interval, the precision 
of the solution cannot be guaranteed but we neverthe- 
less provide the solution over 250 Myr on our Web 
site www. imcce.fr/Equipes/ASD/insola/earth/earth. html 
as reference, and for a possible use, with caution, over 
the full Paleogene period ( up to 65 Myr). 

In order to allow practitioners to test the stability of 
the solution and deduced calibration, we have decided to 
provide the four solutions that have been discussed here, 
La2010a,b,c,d, La2010a being the nominal solution. 

It should be stressed that La2010a is chosen as the 
nominal solution because of its better numerical accuracy, 
but there is no strong clue that La2010a should behave 
better than the other ones. Before —50 Myr, they behave 
practically all in the same way. beyond —50 Myr, the ro- 
bustness of a fit could be tested by changing the solutions. 
Alternatively, in presence of convincing geological record, 
one may conclude that one solution is more probable than 
another one. This will be in some sense a feedback from 
geology to celestial mechanics. 
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Contrary to (Laskar et al. 1993a) and (Laskar et al. 
2004), we have provided here only the eccentricity solu- 
tion. Indeed, although the model for the Earth rotational 
evolution has been improved, the main uncertainty linked 
to the evolution of the tidal dissipative effect in the past 
is still the main unknown parameter for the precession 
and obliquity evolution, and we thus do not believe that 
a new solution would provide more insight than La2004, 
unless a full analysis of the geophysical effects is made, in 
confrontation with the geological records, which becomes 
then out of the scope of the present paper. We thus refer 
to La2004 for precession and obliquity. 

On the other hand, we plan to release soon a full high 
precision, non averaged solution of the rotation and pre- 
cession of the Earth over 1 Myr, as computed from the 
INPOP model. 

After the publication of the La2004 solution (Laskar 
et al. 2004), our goal was to search for an improved so- 
lution, valid over the full cenozoic era. We must say that 
we have not reached this goal. Although the present so- 
lution presents a significant improvement with respect to 
La2004, it is only valid over about 50 Myr. 

The main improvement in the present solution was to 
use a 1 Myr version of the INPOP ephemeris (Fienga et al. 
2008, 2009) in order to fit the initial conditions and pa- 
rameters of the model. As future versions of INPOP will 
appear (Fienga et al. 2010), we will be able to better eval- 
uate the real accuracy of our model. 

At this point, it is still difficult to say wether it will be 
possible to obtain a precise solution over 65 Myr for the 
eccentricity of the Earth. Indeed, in the present solution, 
we have used a more complex model, by adding the five 
main asteroids in the orbital computation, but this added 
also some instabilities in the system, and although we in- 
creased the numerical accuracy of the algorithm by using 
extended precision instead of double precision, we did not 
reach a better numerical accuracy than in La2004. We 
intend to improve on this point, and to increase the nu- 
merical accuracy of the solution, as we would like to be 
sure that the numerical precision is not the limiting factor 
in the final precision of the solution. In order to do so, we 
will also need to improve at the same time on the speed 
of the algorithm, as the present version of our nominal 
solution La2010a took nearly 18 months to complete. 

With the present solution (La2010), we have reached 
the limit of the observational data, and the limit of pre- 
dictability for a precise solution of the orbital evolution 
of the Earth. We have thus decided to provided several 
possible outcomes instead of a single one as usual. The 
solutions La2010a,b,c,d are all available on the Web site 
www.imcce.fr/Equipes/ASD/insola/ earth/ earth, html. 

Practitioners can thus check which of these solutions 
best fit their data beyond —50 Myr. Moreover, it becomes 
clear that the long periodic terms related to 54 — 53 in the 
eccentricity, that appears also as modulation of the the 
405 kyr term in the eccentricity, and the equivalent S4 — S3 
term in the inclination, are some key macroscopic features 
of the orbital solution that are imprinted in the geological 



record. Their precise recovery can thus provide some clue 
for the past chaotic diffusion of the orbital motion of the 
Earth. 
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